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Abstract -We propose a simple microscopic model of molecular dynamics simulation to study 
orientational glass in three dimensions. We present simulation results for mixtures of mildly 
anisotropic particles and spherical impurities. We realize fee solids without orientational order in 
a rotator phase. As the temperature T is lowered, the disordered matrix is gradually replaced by 
four kinds of orientationally ordered, rhombohedral domains. Two-phase coexistence is realized 
in a temperature window. The impurities serve to anchor the orientations of the surrounding 
anisotropic particles, resulting in finely divided domains or medium long-range orientational order. 
We examine the rotational dynamics of the molecular orientations which is slowed down at low T. 
We predict the shape memory effect under a stretching cycle due to inter-variant transformation. 
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Introduction. — Nonsphcrical molecules such as 
KCN can form a crystal without long-range orientational 
order for mild molecular anisotropy ||1 , while liquid crys- 
tal phases can appear for large molecular anisotropy. Such 
crystals are called plastic crystals in a rotator phase. They 
undergo an orientational phase transition as the temper- 
ature T is further lowered, where the crystal structure is 
cubic at high T and non-cubic at low T. With inclusion of 
impurities in such solids, the so-called orientational glass 
has been realized JTj. Around the transitions, a peak in 
the specific heat fTT and softening of the shear modulus 
[l][2] have been observed. In real systems, the molecules 
often have dipolar moments, yielding dielectric anomaly. 
As a similar example, metallic ferroelectric glass, called re- 
laxor, has been studied extensively [3], where frozen polar 
nanodomains have been observed. 

For one-component anisotropic particle systems, there 
have been a number of simulations on the statics [4 ■ 10 



and dynamics TTp2 of the orientational phase transition. 



For two-component anisotropic particle systems, Chong et 
al. examined the slowing-down of the orientational time- 
correlation functions around the glass transition 
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this paper, setting up a simple microscopic model, we will 
investigate the formation of orientational glass. In par- 
ticular, we will examine how impurities can microscopi- 
cally produce orientational disorder, which has remained 
unexplored in the literature. Furthermore, alignments of 



anisotropic particles formig a crystal give rise to lattice 
deformations. As a unique feature in orientational glass, 
heterogeneous strains should emerge on mesoscopic scales. 
In such situations, we may predict soft elasticity and a 
shape memory effect against applied stress. 

Model. — We consider a binary mixture of anisotropic 
particles with number N\ and spherical particles with 
number N%, In this paper, we set N\ + N 2 — 8192. The 
composition of the second species is defined by 



N 2 /(Ni+N 2 ). 



(1) 



We assume relatively small c, so the spherical particles 
may be treated as impurities. The particle positions are 
written as r, (i = 1, ■■■,N). The anisotropic particles 
are assumed to be axisymmetric; then, their orientation 
may be expressed in terms of the solid angles fa and Oi 
(i = l,---,Ni) as rii — (sinflj sin fa, sin#j cos fa, cos Oi). 
The particle sizes of the two species are characterized by 
lengths <T\ and a 2 - The pair potential Uij between par- 
ticles i£a and j G /3 (a,/3 = 1,2) is a truncated mod- 
ified Lennard- Jones potential depending on the particle 
distance r^ = |r,j — rj \ and the directions «.j and rij. For 
r ij > r c — 3o"i it is zero, while for r i? < r c — ?>o~i it reads 



Uu = 4e 



0-12 a 6^ 

(1 + Aij)—^ -g- 



(2) 
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where e is the characteristic interaction energy and 

C Q /3 = (&a + 



^, = 1.1, X=1.2, T = 1 



(3) 



The particle anisotropy is taken into account by the angle 
factor Aij, which depends on the relative direction r,-j = 
r ij ( ri ~ r i) ana - the orientations rii and rij. In this paper, 
we assume the following form, 



X<W> 



(4) 



where x is the anisotropy strength of repulsion. The 8 a \ 
(Spi) is equal to 1 for a = 1 (/3 = 1) and for a = 2 
((3 — 2). In the right hand side, the first (second) term is 
nonvanishing only when i (j) belongs to the first species. 
The Cij ensures the continuity of ?7y at r = r c . 

In our system, the total potential energy is the sum 
U = J2i<i<j<N^ij' while the total kinetic energy K 
arises from the translational velocity fj = dri/dt and the 
rotational velocity hi = drii/dt as 



K = 



E 

Ki<N 



1 



E 5*1* 



(5) 



l<j<JVi 



where mj and mi are the masses, and l\ is the moment 
of inertia of the first species. We set m\ — mi in our 
simulation. The molecular rotation around the symmetry 
axis parallel to rii does not change U, so we may neglect 
its kinetic energy. The Newton equations of motion for 
translation and rotation are written as 1141 



dU 



or 

hrii x h. t = -rii x 



(i = l,---,N), 
dU 



(6) 



(i = l,---,N 1 ), (7) 



where 7% = d 2 ri/dt 2 and hi = d 2 rii/dt 2 . Hereafter, we 
measure space, time, and temperature in units of <T\, tq = 
o\ ^mi/e, and e/fcs, respectively. 

Treating equilibrium or at least nearly steady states, 
we furthermore attach a Nose-Hoover thermostat 
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to 

all the particles. That is, we added the thermostat terms 
— C^a(t) m a r i and — Cmi(t)IiniXfii in the right hand sides 
of Eqs.(6) and (7), respectively, where Cnh(0 obeys 



Tnh^tCnhW 



■l + KftkBTfiNfi + N!)], (8) 



with tnh being a short relaxation time taken to be 0.1. 

We may envisage the anisotropic particles as spheroids 
depending on the parameter \ from Eq. (4) . For particles i 
and j of the first species, minimization of Uij with respect 
to r t j yields r V] = 2 1 / 6 (l + A i j) 1 / 6 (r 1 . For \ > 0, the short 
and long diameters, a s and ag, are given by 



2 1/6 (Ti, a, = (l + 2 X ) 1 /6 2 V< 



(9) 



for the perpendicular and parallel orientations of rii and 
rij with respect to fy, respectively. The inertia momen- 
tum of the first species in Eq.(7) is given by 



h = (of + a 2 s )mi/20. 



(10) 




c=0.3 



Fig. 1: Frozen domain structures composed of four rhombohe- 
dral variants for 02/01 = 1.1, x — 1-2, and T = 0.1 at fixed 
volume. With increasing the composition c of the impurities 
(black points), the orientational disorder increases, leading to 
a decrease in the domain size. The particle color is blue, green, 
and red for rii = (±1, 0, 0), (0, ±1, 0), (0, 0, ±1), respectively. 



In our simulation, we set <Jijo\ = 1.1 and x = 1-2. 
From Eq.(9) the aspect ratio is ai/a s — 1.23. For these 
mild parameter values, we realized fee plastic solids with- 
out long-range orientational order at relatively high T. 
For large Xi sa Y 10> we realized liquid crystal phases in 
our model. It is worth noting that our angle-dependent 
potential (2) is analogous to the Gay-Berne potential for 
anisotropic molecules used to simulate mesophases of liq- 
uid crystals [16] . Similar angle-dependent potentials have 
also been used for lipids forming membranes. 
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Four-variant states at fixed volume. Starting 
with liquid states at T — 2, we quenched the system to 
T = 0.5 below the melting and waited for 10 to real- 
ize fee plastic crystals. Here, a few stacking faults were 
formed in most runs. We then cooled the system to a 
final temperature and waited for 3 x 10 4 , where orienta- 
tional order developed. In Figs. 1-7, we fixed the system 
volume and shape imposing the periodic boundary condi- 
tion. In terms of the molecular volumes V\ — ira 2 ag/6 and 
V2 = 7t2 1 / 2 (t|/6, the packing fraction <^ pac k is given by 



(N lVl + N 2 v 2 )/V. 



(11) 



We set </> pac k = 0.82 in Figs. 1-7. Then the system length 
is about 20 and the pressure is about ie/af. 

In Fig.l, the orientational domains can be seen at 
T = 0.1, where the particle positions and orientations are 
nearly frozen in time. For c = 0, the system undergoes 
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Simulation of orientational glass 




Fig. 2: Anchoring of spheroidal particles around a single impu- 
rity (left) and two associated impurities (right) for a%lo\ = 1.1, 
X = 1.2, and T = 0.1. 



a structural phase transition, where the lattice structure 
changes from fee to rhombohedral almost without dilation 
change. The lattice constant is equal to a = 1.65 through 
the transition. The ordered phase consists of four rhom- 
bohedral variants where the angles of the rhombuses are 
83° and 97° . The separation distance of the closely packed 
(111) planes is increased from <z/y3 = 0.952 to 1.08, be- 
cause the orientation vectors rti are aligned in the [111] 
direction. With increasing c, the orientational disorder is 
gradually increased, leading to pinning of finely divided 
domains. We eventually obtain a glass state at c — 0.3. 

In Fig. 2, the anisotropic particles are aligned in the 
perpendicular directions (_L f^) around one or two im- 
purities. This parallel anchoring disturbs the orientation 
order. In Fig. 3, in a (111) plane, we display the four rhom- 
bohedral variants with different orientations for c = 0.2. 
The interfaces between the variants tend to be trapped at 
clustered impurities. We also notice that the junction an- 
gles, at which two or more domain boundaries intersect, 
are approximately multiples of 7r/6. Similar patterns were 
observed on hexagonal planes in a number of experiments 
on alloys after structural phase transitions 
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Impurity clustering. — Figure 4 shows significant 
impurity clustering, which appeared during quenching 
from liquid. The average impurity number per cluster is 
5.5 for c = 0.1 and there is a big percolated cluster for 
c = 0.2. Here, two impurities are treated to belong to the 
same cluster if their distance is smaller than 1.3. In the 
snapshot of c = 0.3 in Fig.l, the impurity clustering is 
closely related to the heterogeneity in the orientations. 

To explain this effect, we compare the solvation en- 
ergy of a single impurity and that of two associated im- 
purities (see Fig. 2). In terms of the potential energy 
Uj = Ujk/2 of particle j, it is estimated as 



sol 



E (Uj-Vi)- 



(12) 



nearby j£l 



Here, £/i mp is the contribution from the impurities un- 
der consideration, the summation is over the nearby 
anisotropic particles j with ry < 1.3 (i 6 2 and j € 1), 
and U\ is the average potential energy of the non-neighbor 
anisotropic particles separated from any impurities longer 
than 1.3. At T = 0.1, U so \ is calculated as = — 6.3e 



c=0.2 a 2 /cr 1 =11,X=1.2, T=0.1 




(111) plane 
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Fig. 3: Crosssectional snapshot in a (111) plane, where four 
rhombohedral domains composed of spheroidal particles can be 
seen with impurities (black) disturbing the orientations. Here, 
c = 0.2, <j 2 /o-i = 1.1, x = 1-2, and T = 0.1 . 



for a single impurity and U^) = — 13. 8e for a dimer, where 



sol 

-4.5e. The difference AJ7 sol = U<$- 



2C7«(~-1.2e) 
is the association energy. At T — 0.25, we have AU so \ ~ 
— 1.8e. The total potential energy is thus lowered with 
impurity clustering in the present case. 

Coarse-grained orientation order parameter. 

For each particle i of the first species, we introduce the 

orientation tensor Q i = {Qi^ u } (fJ,, v = x, y, z) as 



1 



UiTl, 



E ' 

j'G bonded 



1 <-> 

-3 1 ' 



(13) 



where 1 = {5^} is the unit tensor. The summation is 
over the bonded particles of the first species with |ry| < 
1.5o"i and is the number of these particles. If a fee 
lattice is formed, it includes the nearest neighbor particles, 
so n l h ~ 12. We define the amplitude of the orientational 
order for each anisotropic particle i as 



2 



(14) 



Here, Si ~ 0.1 in disordered regions due to the thermal 
fluctuations, but it increases up to unity within rhombo- 
hedral domains at low T. In Fig. 5, the average (S) = 
J2i<i<N! Si/Ni represents the overall degree of orienta- 
tion order. In our simulation, we realized only uniaxial 

states, where we have Qi= S i (didi— 1 /3) in terms of 
the amplitude Si and the director d t . 



Coexistence of high and low temperature phases. 

— In Fig. 5, the degree of orientational order (S) in- 
creases continuously as T is lowered. For small c at fixed 
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1.2, T=0.1 




Fig. 4: Snapshots of impurities composing clusters with sizes 
> 5 for c = 0.1 (left) and for c = 0.2 (right) with mesoscopic 
heterogeneities. The data are the same as those in Fig.l. 
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Fig. 5: (S) = J2 ie i S i/ N * vs T > for c = °> 0- 1 ' °- 2 > °- 3 > and °- 4 
for a%ja\ = 1.1, and \ — 1.2. This quantity represents the 
overall degree of orientational order. 



c=o.2, o 2 /G 1 =^.^,% = ^.2 




T=0.15 



Fig. 6: Left: Distribution of the orientation amplitude P(S) in 
Eq.(15) for various T, where c = 0.2, 02/01 = 1.1, and \ = 1.2. 
It exhibits double peals for T = 0.13, 0.15, and 0.18 repre- 
senting coexistence of cubic and rhombohedral regions. Right: 
Snapshot of anisotropic particles with Si > 0.6 at T = 0.15, 
forming ordered domains embedded in a disordered matrix. 



volume and shape, however, we find coexistence of cu- 
bic and rhombohedral regions in a temperature window, 
which is roughly given by 0.26 <T <0.29 for c = and 
0.14 <T <D.19 for c = 0.2. To examine this coexistence, we 
introduce the distribution of the orientational amplitude, 



P(S) 



1 



E 

l<i<Nl 



Si)}- 



(15) 



Taking the average (• • •) over six runs we obtained P(S) 
in Fig. 6, which exhibits two peaks for 0.13 < T < 0.18 at 
c = 0.2. We also give a snapshot of the ordered anisotropic 
particles with Si > 0.6. For c = 0, the interfaces between 
the ordered and disordered regions are temporally fluctu- 
ating with small amplitudes. Small amounts of impurities 
can pin the interfaces and the ordered regions are stabi- 
lized far from the impurity clusters. 

Rotational dynamics. — In crystal and glass, elon- 
gated particles often undergo turnover motions, rii — > 
— rij. These events take place in a microscopic time 
(~ 1), while the characteristic time t\ between succes- 
sive turnover motions grows at low T. Let us consider the 
following angle relaxation function, 



1 

Ad 



E 

l<i<ATi 



5(C-ni(t + to)-ni(t ))), (16) 



where the average (• • •) is taken over the initial time to 
and over six runs in this paper. Furthermore, in terms 
of the ^-th order Legendre polynomials Pi, we may define 
the £-th order rotational relaxation functions as 
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c t {t)= / dCP e (C)W((,t). 



(17) 



Here C\(t) decays on the timescale of n, while Ci{i) is 
unchanged by the turnover motions. In orientational glass, 
the ultimate relaxation of C2 (t) is due to the orientational 
configuration changes involving the surrounding particles, 
so its relaxation time T2 much exceeds t\ . 

In Fig.7, we plot C x {t) and C 2 (t) vs t for c = 0.2. For 
t <1 they relax considerably due to the thermal motions. 
The decay of C\{t) is slower for lower T, while C2(t) appar- 
ently tends to a positive constant f% for large t, as in the 
previous simulation [13] . In our case, this plateau appears 
because the anchoring of the anisotropic particles around 
the impurities becomes nearly permanent at low T. Thus 
fi increases with lowering T. In Fig.7, we also present 
the time-evolution of the relaxation function W((,t) for 
t = 200, 2000, and 20000, which is peaked at C = ±1- The 
peak height at f = — 1 increases in time from 0. Also 
shown is the turnover probability defined by 



Wto(i) 



-l+AC 



d(W(c,t), 



(18) 



where we set A£ = 0.2. We notice that W to (t) grows 
linearly in time in the very early stage. Thus we may 
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Fig. 7: Orientation relaxation functions C\ (t) in (a) and C2 (t) 
in (b) for c = 0.2, ai/ai = 1.1, and x = 1.2. Ci(t) decays 
due to turnover motions and its relaxation time t\ is marked 
on each curve (see Eq.(19)). C%{£) can relax only due to the 
orientational configuration changes, (c) Distribution W(£, t) 
in Eq.(16) and (d) turnover fraction in Eq,(18) at three times. 



define Ti = T\ (T) from the following linear form, 

Wt (i) = rf 1 *, (19) 

which holds for i/n <D.l. For 0.1 <£/ti <1, the probabil- 
ity of multiple turnovers becomes appreciable resulting in 
a deviation from the linear behavior (19). However, the 
scaling form W to (t) = /to(^/ r i) fairly holds for t/n < 1. 
For i/n >4, the orientational structural relaxation comes 
into play and the scaling in terms of t\ docs not hold. 

Shape memory effect. — In the presence of multi- 
variant orientational order, the shape-memory effect (me- 
chanical hysteresis) emerges almost without dislocation 
formation. This effect is well-known for shape memory 
alloys such as Ti-Ni [20], where a structural phase tran- 
sition is caused by the relative atomic displacements in 
each unit cell. In our case, under stretching along the z 
axis, domains with rii nearly parallel to the z axis grow 
yielding an increase in the strain. We are not aware of any 
experiments on this effect in anisotropic particle systems. 



In this mechanical effect, the system shape changes 
Thus we assumed a Parrinello- Rahman barostat MpT 



as 

well as the Nose-Hoover thermostat. Before streching we 
prepared a multi-variant initial state at t = 0, where the 
packing fraction </> pac k m Eq.(ll) was 0.75 and the pressure 
was zero. For t > wc applied the stress a a = — (H zz ), 
setting (Hxx) = (n w ) = 0, where 11^ are the stress com- 
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A)a a =0.35 



B) 0.45 



Fig. 8: Shape memory effect under stretching along the z axis 
in orientational glass with c = 0.3 in stress-controlled simula- 
tion, where T = 0.05, di/cri = 1.1, and \ = 1-2. Top left: 
Strain e z vs stress a a , where e z increases slowly in the hard 
ranges < a a < 0.38 and 0.42 < a a < 0.5 and steeply in the 
soft range 0.38 < a a < 0.42. The return path is reversible with 
large Young's modulus. After this cycle, the residual strain is 
0.1, which vanishes upon heating to T = 0.2. Top right: Time- 
evolution of four variant fractions on the stretching path. For 
o a > 0.38 the two variants with their orientations rii nearly 
parallel to the z axis become dominant. Bottom: Snapshots of 
the particle configurations at (A) a a = 0.35 and (B) a a = 0.45. 



ponents with (• • •) being the space average. We allowed 
the system to take a rectangular shape. In terms of the 
system length L z (t) along the z axis, the average strain is 



= L z (t)/L z (0)-l. 



(20) 



We define effective Young's modulus E e by de z /d<r a = 
1/E e . It follows the effective shear modulus /x e from the 
formula E e = 3fi e / (l+/j, e /3K) in classical elasticity, where 
K is the bulk modulus. In our case K is about 50 and is 
much larger than fj, e , so E e = 3/i e . Hereafter, a a and E e 
will be measured in units of e/af. 

In Fig. 8, we first increased a a at da a /dt = 5 x 10~ a 
from up to 0.5 at T = 0.05. In the left panel, we find 
E e - 17 for < a a < 0.38, E e ~ 3 for 0, 38 < cr Q < 0.42, 
and E e ~ 20 for 0,42 < a a < 0.5. In the second range, 
the solid is soft because the variants elongated along the 
z axis grow due to inter- variant transformation. Here, the 
orientation vectors rii — (n^, ni y , m z ) of the four variants 
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are roughly given by (a) (0.3,-0.1,0.9), (b) (-0.7,0.1,0.7), 
(c) (0.7,0.7,0.2), and (d) (0.5,-0.9,0.1) in Fig.8. In the 
right panel, the variant (a) grows, the variant (b) remains 
nonvanishing, and the variants (c) and (d) disappear. Sec- 
ondly, we decreased a a back to at da a /dt = —5 x 10~ 5 . 
In this return path, E e was about 20, the disfavored vari- 
ants no more appeared, and a remnant strain of order 0.1 
remained. Thirdly, at a a — 0, we increased T from 0.05 
to 0.2 into the disordered phase. After this heating, the 
system became isotropic. 

Summary. — We have siudied the orientational glass 
using the potential (2). Due to the anisotropic factor A^, 
the anisotropic particles have the aspect ratio (1 + 2x) 1 ^ 6 - 
In our simulation, the anisotropy parameter x is 1.2 and 
the size ratio ct^/cfi is 1.1. As a result, fee plastic crystals 
appear at relatively high T, which undergo the orientation 
transition at lower T. The impurities serve to pin the 
orientations of the surrounding anisotropic particles. 

We briefly summarize our main results. In Fig.l, we 
have illustrated four-variant ordered states for various c 
at T = 0.1. With increasing c, rhombohedral domains are 
finely divided and their typical size is decreased. In Fig. 4, 
we have shown that the impurities are heterogeneously 
distributed. In Fig. 5, we have plotted the average orien- 
tational amplitude (S), which increases continuously with 
lowering T. In Fig. 6, we have shown large-scale coexis- 
tence of the disordered and ordered phases in a tempera- 
ture window. In Fig. 7, the orientational time-correlation 
functions have been plotted for c = 0.2, where C\{t) de- 
cays due to turnover motions. In Fig.8, we have examined 
the shape memory effect in orientational glass. 

We next make some remarks below. 

(i) For small c, we need to know how the characteristic do- 
main size is determined. For moderate impurity concen- 
trations (c >0.2 in this work), the orientational disorder is 
proliferated. These aspects should be further studied. 

(ii) The phase transition delicately depends on whether 
simulations are performed at fixed volume or at fixed 
stress. Some salient features at fixed stress are as follows, 
(a) For c = 0, a first-order phase transition occurs with a 
shape change, (b) For c > 0, both multi-domain states and 
single-domain states can be realized for the same param- 
eters as in Fig.8. (c) For c > 0, the two-phase coexistence 
in Fig. 6 can still be realized even at fixed stress. 

(iii) In our simulation, the system keeps the crystalline or- 
der. However, for larger size ratio (T2/&1, sa y 1-4, we found 
an increase in the positional disorder leading to polycrystal 
and positional glass. Competition of positional and orien- 
tational glass transitions can then be studied. For large 
anisotropic parameter x, we may also study complicated 
phase behavior of two-component liquid crystals. 

(iv) In this paper, the impurities are spherical and slightly 
larger than the anisotropic particles. By modifying the 
potential form, we may also treat other types of impu- 
rities. For example, their attractive interaction with the 
host anisotropic particles may be anisotropic; then, the 



anchoring can be homeotropic. 
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